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Abstract 

Energies and spatial observables for the 1 E g state of the hydrogen molecule in mag- 
netic fields parallel to the proton-proton axis are calculated with a guided random walk 
Feynman-Kac algorithm. We demonstrate that the accuracy of the results and simplicity 
of the method may prove it a viable alternative to large basis set expansions for small 
molecules in applied fields. 

Suggested PACS numbers: 31.15 +q, 31.15 -Kb, 33.15 Bh 

1. Introduction 

The physics of the H2 molecule in magnetic fields is a problem of interest in molecular 
physics [1,2], astrophysics [3], and to some extent in condensed matter physics [4]. A series 
of papers examining Hi in magnetic fields parallel to the proton-proton axis have recently 
been published [5,6]. These papers use basis set expansions to compute potential energy 
curves and equilibrium values of the proton-proton radial separation for several of the low 
lying states of the molecule. While this program has been succesful, it is useful to have 
alternate methods that have their own advantages. Work has been done on this problem 
with a fixed phase Monte Carlo algorithm [7] that is general and may prove to be the 
eventual method of choice for other molecules in magnetic fields. Here however we have 
chosen a guided random walk Feynman-Kac (GRWFK) algorithm [8,9,10] because of it's 
ease of implementation and extension to other molecular systems. GRWFK also has the 
advantage of avoiding trial wavefunction contamination of quantities calculated with the 
method. 

This brief report is organized as follows: In section 2 we give an outline of the calcu- 



lational method. In section 3 we apply the method to the 1 T, g state of H 2 in a magnetic 
field. Section 4 presents results for the energies and electronic (Vr*) values of the molecule 
for some representative fields. Section 5 is reserved for conclusions. 

2. Calculational method 

We give a brief discussion of the calculational method for one dimensional systems. 
The extension to higher dimensionality is very simple and is an advantage of the formalism. 
The Feynman-Kac formula [11,12] for a walker starting a walk at x', T' = 0, is 
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U(x, T, x 1 = 0, T' = 0) = (exp[ - I ds V(B(s) +*)]> = £ ^ n (x)^ n (0) e - E - T . (1) 
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U(x,T, 0,0) is the Euclidian quantum mechanical propagator and the brackets denote an 
expectation value over the Weiner measure [13]. The argument of the exponential is an 
integral of the potential along a Gaussian process of mean zero and variance one [14] (we 
choose this as a Bernoulli sequence of plus and minus ones) . Computationally this amounts 
to evaluating 
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a large number N times. In practice the energy is better approximated by 
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with Ti,T 2 , two sufficiently large times. 

It is known the F-K formula as written in eq.(l) and it's discretized form in eq.(2) are 
formally correct but plagued by slow convergence and large variance when employed for 
numerical work. There are many schemes in the literature that address this problem. Here 
we adopt a method that has the advantage of simplicity and adaptability to parallelization. 

An unconstrained random walk has equal probability P = | to step left or right. The 
probability distribution p(x, r + At) for the arrival of a walker at x, r is 

p(x, t + At) = ^p(x - Ax, t) + + Ax, r). 



This leads to 

dp 1 h 2 d 2 p 
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dr 2 h T d 2 x' 

If instead the walk probabilities for a step right or left are modified by a guiding function 
g(x) such that 

, 2 

then in the continuum limit with -ff > 1 a modified diffusion equation (inclusion of the 

potential V is simple [8]) for the probability distribution becomes 

dp ld 2 p dg dp d 2 g 
-d-r=—2W- X - ^xWx + P Wx + Vp - (6) 

Although formally p ought to be replaced with another symbol to account for inclusion of 
V(x) it should cause no confusion here. 
By setting p = e~ 2g U eq.(6) becomes 

dU ld 2 U „ TT 

n( x) = v( x )-2\ d ^] + d -^l. (8) 
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Since g = —\lnU (from p oc \1/ 2 ) knowledge of the solution of eq.(6) give would the 
exact guiding function g(x) and vice- versa. Of course U is not generally known, but by 
choosing g(x) to incorporate what is known about the general character of the solution the 
variance of calculated observables can be reduced and convergence to the final result can be 
substantially improved. In what follows we evaluate eq.(7) with conditional probabilities 
given by eq.(5), and use Q(x) rather than V(x) as the argument of the exponential in eq.(2). 
For the ground state U oc *o so we refer to trial wavefunctions rather than propagators 
below. 

The ground state expectation value of an operator 0(ti) can be found by evaluating 

{0 ,_ E ri O(ri)|*(ri)| 2 
E ri l^(ri)| 2 

For large times Ti, T 2 the rightmost term of eq. (1) can be used to generate |\P(ri)| , hence 
(O): Evolve the walk to to time Ti, tallying the value of the exponential of eq.(2) for each 



path at Ti. Multipy this value by (0(Ti)),then evolve the walk to T 2 . With 

T 2 
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this yields 



<0> 



It is worth noting that we do not generate random walks on a three dimensional grid for 
each particle but rather as three separate one dimensional grids. This is why the three 
dimensional spatial branching factor of 6 used by Barnes et.al. does not appear above. 

3. Application to H 2 in a magnetic field. 

In the clamped nuclei approximation the Hamiltonian of the if 2 molecule in a magnetic 
field B may be written in natural units as 
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with £? = 2.3505 x 10 9 G. Taking B along the z axis reduces eq. (9) for the 1 E £/ state to 
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For the trial 1 E 5 wavefunction it is convenient to first define the auxilary quantities 

f* = -\*±*\- 



The trial wavefunction is then 
* t (ri,r 2 ,R) = 



oh + _|_ e /r 



l(xi 2 +m 2 ) 



x 
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(11) 



It is straightforward to insert the of eq. (11) into the expression for g to obtain 
f2 and the modified walk probabilities. We note in passing that this choice amounts to a 



baseline level trial function. No a priori attempt was made to include free parameters for 
later optimization. 

4. Results 

In table 1 we give E in atomic units for several values of 7 at the indicated R eq . We 
also give Eq for 7 = 2.127207 at R eq = 1.07 a.u. There are two values in the literature for 
Eq at this field strength [6,15]. Our value confirms the result of Detmer et. al. in [6]. 

These results show good agreement with [6] for the relatively small number of sample 
paths considered here. The value n = 800 was arrived at by starting at n = 200 and 
increasing n by 100 steps per run for a few thousand sample paths, n = 500 was found to 
be sufficient for all 7 considered. 300 more steps were added as a check on the stability 
of convergence. We give the results for the larger n. A similar method was employed for 
the determination of Ti and T 2 . (We did not include disassociation values since they were 
also the same as those given in [6].) In table 2 we show results for calculations of (Vr^) 
for some values of 7. We were interested in the value of the field for which the rms value 
equaled the proton's separation believing this to be a reasonable albeit naive measure of 
the onset of the transition to the free atomic limit. Clearly the field value at which this 
would occur is large enough to align the spins making the triplet rather than the singlet 
state the correct state to investigate. It should be noted that the times Ti and T 2 that 
give convergence of E usually were not sufficient to converge (Vr^). This is reflected in 
the larger values of T\ and Ti in table 2. 

5. Conclusions 

The accuracy of the GRWFK method for this problem suggests that faster convergence 
with basis set expansions might be achieved with LCAO trial wavefunctions functions 
multiplied by oscillator eigenstates. This is a direction best explored by those working 
already with that method. Certainly basis sets expansions prove (at least for fixed nuclei) 
a superior method for digit aquisition than Monte Carlo methods. The trade off between 
accuracy and ease of implementaion for molecules like H3 or even Hs + in a magnetic field 
is an open question. 

Calculating the properties of low-lying excited states is not unfeasable but is known 



to be limited by the (lack of) knowledge of the nodal structure of the state [16,17]. The 
fixed node method [18,19] has been shown as a good approximation for small molecular 
systems [20]. For biased random walks the nodal structure can be built into the the trial 
functions although not without increased complexity. Clearly states of non-zero orbital 
angular momentum will have to be considered by incorporating the techniques discussed 
in [7] to GRWFK. This is a topic of current investigation. 
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Table 1 1 E £/ state energies as functions of magnetic field 
B and R eq . All units are in a.u. Bracketed numbers are 
last digit uncertainty estimates from variance calculations 
for N = 50000 sample paths, n = 800, T x = 7, T 2 = 8 for 
all entries. 
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Table 2 1 E 3 \fr* values in a.u. for some values of B. Data 
below are for n = 800, T x = 10, T 2 = 12 and N = 50000 
sample paths. Larger times are needed to converge Vr* 
compared to those needed for Eq. The last entry in the 
table serves as a check on the dissacociative limit. 





Req 




.01 


1.40 


1.565(3) 


1.0 


1.24 


1.326(2) 


2.0 


1.09 


1.141(3) 


5.0 


0.86 


0.882(2) 


7.5 


0.74 


0.782(1) 


10.0 


0.70 


0.658(2) 


0.0 


50.0 


1.495(2) 



